Measurement exploratory data analysis¶
DEVICES = [
'beaglebone-fan',
'beaglebone-compressor',
'beaglebone-pump',
'beaglebone-refrigerator',
'mafaulda-a',
'mafaulda-b'
]
DEVICE = DEVICES[0]
T_WAVEFORM = 5 # (1 = MaufaulDa, 5 = others)
T_SEC = T_WAVEFORM
NFFT = 2**10 # (2**10 = MaufaulDa, 2**14 = others)
F_LIMIT = None # (3000 = MaufaulDa, None = others)
import os
import pandas as pd
import numpy as np
import matplotlib.pylab as plt
from tabulate import tabulate
from IPython.display import Markdown, HTML
from tqdm.notebook import tqdm
from typing import List, Tuple
from scipy.signal import find_peaks, butter, lfilter
from tsfel.feature_extraction.features import fundamental_frequency
from zipfile import ZipFile
import sys
sys.path.append('../')
from vibrodiagnostics import mafaulda, selection, discovery, models
def beaglebone_measurement(filename: str, fs: int) -> Tuple[str, pd.DataFrame]:
g = 9.81
milivolts = 1800
resolution = 2**12
columns = ['x', 'y', 'z']
ts = pd.read_csv(filename, delimiter='\t', index_col=False, header=None, names=columns)
# Calculate amplitude in m/s^2 Beaglebone Black ADC and ADXL335 resolution (VIN 1.8V, 12bits)
for dim in columns:
ts[dim] = ts[dim] * (milivolts / resolution) # ADC to mV
ts[dim] = (ts[dim] / 180) * g # mV to m/s^2 (180 mV/g)
ts[dim] -= ts[dim].mean()
ts['t'] = ts.index * (1 / fs)
ts.set_index('t', inplace=True)
return (os.path.basename(filename), ts, fs, ts.columns) # last is feature columns
def beaglebone_dataset(filenames: List[str], fs: int) -> List[Tuple[str, pd.DataFrame]]:
dataset = []
for filename in filenames:
name, ts, fs, cols = beaglebone_measurement(filename, fs)
dataset.append((name, ts))
return dataset
def lowpass_filter(data, cutpoint, fs, order=5):
b, a = butter(order, cutpoint, fs=fs, btype='lowpass')
y = lfilter(b, a, data)
return y
def mafaulda_dataset(
place=['ax', 'ay', 'az'],
features_path = '../../datasets/features_data/',
mafaulda_path='../../datasets/MAFAULDA.zip',
rpm=2500,
lowpass_hz=10000):
metadata_filename = os.path.join(features_path, selection.MAFAULDA_METADATA)
faults = {
'normal': 'normal',
'imbalance': 'imbalance',
'horizontal-misalignment': 'misalignment',
'vertical-misalignment': 'misalignment',
'overhang-cage_fault': 'cage fault',
'underhang-cage_fault': 'cage fault',
'underhang-ball_fault': 'ball fault',
'overhang-ball_fault': 'ball fault',
'overhang-outer_race': 'outer race fault',
'underhang-outer_race': 'outer race fault',
}
metadata = pd.read_csv(metadata_filename, index_col='filename')
metadata.reset_index(inplace=True)
metadata = models.fault_labeling(metadata, faults)
files = pd.DataFrame()
# Worst severity and mid rpm
for name, group in metadata[metadata['rpm'] >= rpm].groupby(by='fault', observed=False):
files = pd.concat([
files,
group[
group['severity_level'] == group['severity_level'].max()
].sort_values(by='rpm', ascending=True).head(1)
])
ordering = {
'normal': 0,
'misalignment': 1,
'imbalance': 2,
'cage fault': 3,
'ball fault': 4,
'outer race fault': 5,
}
source = ZipFile(mafaulda_path)
dataset = len(files) * [0]
for index, file in files.iterrows():
ts = mafaulda.csv_import(source, file['filename'])
ts = ts[place]
ts.columns = ts.columns.str.extract(r'(\w)$')[0]
for axis in ts.columns:
ts[axis] = lowpass_filter(ts[axis], lowpass_hz, file['fs'])
pos = ordering[file['fault']]
dataset[pos] = ((file['fault'] + ' (' + file['filename'] +')', ts))
return dataset
Load dataset
if DEVICE == 'beaglebone-fan':
Fs = 2500
path = '../../inspections/fan/'
files = [
'1_still.tsv', '2_still.tsv', '3_still.tsv',
'1_up.tsv', '2_up.tsv', '3_up.tsv',
'1_down.tsv', '2_down.tsv', '3_down.tsv'
]
files = [os.path.join(path, name) for name in files]
DATASET = beaglebone_dataset(files, Fs)
elif DEVICE == 'beaglebone-compressor':
Fs = 2500
path = '../../inspections/datacentres/shc3/'
files = [
'k3_1.tsv', 'k3_2.tsv', 'k3_3.tsv', 'k3_4.tsv',
'k5_1.tsv', 'k5_2.tsv', 'k5_3.tsv', 'k5_4.tsv'
]
files = [os.path.join(path, name) for name in files]
DATASET = beaglebone_dataset(files, Fs)
elif DEVICE == 'beaglebone-pump':
Fs = 2500
path = '../../inspections/bvs/'
files = [
'bvs_1_hore.tsv', 'bvs_2_hore.tsv'
#, 'bvs_3_motor.tsv', 'bvs_4_motor.tsv'
]
files = [os.path.join(path, name) for name in files]
DATASET = beaglebone_dataset(files, Fs)
elif DEVICE == 'beaglebone-refrigerator':
Fs = 2500
path = '../../inspections/home-refrigerator/'
files = [
'ch1.tsv', 'ch2.tsv', 'ch3.tsv', 'ch4.tsv', 'ch5.tsv'
]
files = [os.path.join(path, name) for name in files]
DATASET = beaglebone_dataset(files, Fs)
elif DEVICE == 'mafaulda-a':
Fs = mafaulda.FS_HZ
DATASET = mafaulda_dataset(place=['ax', 'ay', 'az'])
elif DEVICE == 'mafaulda-b':
Fs = mafaulda.FS_HZ
DATASET = mafaulda_dataset(place=['bx', 'by', 'bz'])
for name, ts in DATASET:
display(Markdown(f'**{name}**'))
ts.info()
print()
1_still.tsv
<class 'pandas.core.frame.DataFrame'> Index: 38200 entries, 0.0 to 15.2796 Data columns (total 3 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 x 38200 non-null float64 1 y 38200 non-null float64 2 z 38200 non-null float64 dtypes: float64(3) memory usage: 1.2 MB
2_still.tsv
<class 'pandas.core.frame.DataFrame'> Index: 38200 entries, 0.0 to 15.2796 Data columns (total 3 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 x 38200 non-null float64 1 y 38200 non-null float64 2 z 38200 non-null float64 dtypes: float64(3) memory usage: 1.2 MB
3_still.tsv
<class 'pandas.core.frame.DataFrame'> Index: 38200 entries, 0.0 to 15.2796 Data columns (total 3 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 x 38200 non-null float64 1 y 38200 non-null float64 2 z 38200 non-null float64 dtypes: float64(3) memory usage: 1.2 MB
1_up.tsv
<class 'pandas.core.frame.DataFrame'> Index: 38200 entries, 0.0 to 15.2796 Data columns (total 3 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 x 38200 non-null float64 1 y 38200 non-null float64 2 z 38200 non-null float64 dtypes: float64(3) memory usage: 1.2 MB
2_up.tsv
<class 'pandas.core.frame.DataFrame'> Index: 38200 entries, 0.0 to 15.2796 Data columns (total 3 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 x 38200 non-null float64 1 y 38200 non-null float64 2 z 38200 non-null float64 dtypes: float64(3) memory usage: 1.2 MB
3_up.tsv
<class 'pandas.core.frame.DataFrame'> Index: 38200 entries, 0.0 to 15.2796 Data columns (total 3 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 x 38200 non-null float64 1 y 38200 non-null float64 2 z 38200 non-null float64 dtypes: float64(3) memory usage: 1.2 MB
1_down.tsv
<class 'pandas.core.frame.DataFrame'> Index: 38200 entries, 0.0 to 15.2796 Data columns (total 3 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 x 38200 non-null float64 1 y 38200 non-null float64 2 z 38200 non-null float64 dtypes: float64(3) memory usage: 1.2 MB
2_down.tsv
<class 'pandas.core.frame.DataFrame'> Index: 38200 entries, 0.0 to 15.2796 Data columns (total 3 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 x 38200 non-null float64 1 y 38200 non-null float64 2 z 38200 non-null float64 dtypes: float64(3) memory usage: 1.2 MB
3_down.tsv
<class 'pandas.core.frame.DataFrame'> Index: 38200 entries, 0.0 to 15.2796 Data columns (total 3 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 x 38200 non-null float64 1 y 38200 non-null float64 2 z 38200 non-null float64 dtypes: float64(3) memory usage: 1.2 MB
for name, ts in DATASET:
display(Markdown(f'**{name}**'))
display(tabulate(ts.describe(), headers='keys', tablefmt='html'))
ts.boxplot(grid=True)
plt.show()
1_still.tsv
| x | y | z | |
|---|---|---|---|
| count | 38200 | 38200 | 38200 |
| mean | 4.53297e-15 | 3.95579e-15 | 3.23185e-15 |
| std | 0.874418 | 1.2639 | 0.0920418 |
| min | -1.83871 | -2.52106 | -0.348954 |
| 25% | -0.5933 | -1.108 | -0.061552 |
| 50% | 0.00545463 | -0.00628912 | 0.0102986 |
| 75% | 0.580259 | 1.11937 | 0.058199 |
| max | 1.96937 | 2.43663 | 0.345601 |
2_still.tsv
| x | y | z | |
|---|---|---|---|
| count | 38200 | 38200 | 38200 |
| mean | 8.66044e-16 | 6.06082e-15 | 2.94447e-15 |
| std | 2.56509 | 2.48179 | 0.322518 |
| min | -4.5757 | -3.73282 | -0.952338 |
| 25% | -2.44414 | -2.63111 | -0.281733 |
| 50% | 0.429887 | 0.0752644 | 0.00566968 |
| 75% | 2.27405 | 2.56608 | 0.269122 |
| max | 4.04637 | 3.66779 | 0.987628 |
3_still.tsv
| x | y | z | |
|---|---|---|---|
| count | 38200 | 38200 | 38200 |
| mean | -8.10354e-15 | -3.27184e-16 | 1.05735e-14 |
| std | 2.15156 | 2.64354 | 0.396043 |
| min | -3.02308 | -4.12174 | -0.786905 |
| 25% | -2.13692 | -2.58893 | -0.379752 |
| 50% | -0.244853 | -0.0502045 | 0.00345146 |
| 75% | 1.95857 | 2.56635 | 0.362704 |
| max | 3.92248 | 4.21293 | 0.841708 |
1_up.tsv
| x | y | z | |
|---|---|---|---|
| count | 38200 | 38200 | 38200 |
| mean | -2.19487e-17 | 1.16068e-15 | -4.63229e-15 |
| std | 0.850544 | 1.19633 | 0.160765 |
| min | -2.0824 | -2.62078 | -0.818255 |
| 25% | -0.381937 | -0.704763 | -0.0757992 |
| 50% | 0.00126648 | -0.010207 | -0.00394865 |
| 75% | 0.40842 | 0.756199 | 0.0679019 |
| max | 2.01308 | 2.57641 | 0.930109 |
2_up.tsv
| x | y | z | |
|---|---|---|---|
| count | 38200 | 38200 | 38200 |
| mean | -3.61707e-15 | -5.37892e-15 | 3.77964e-16 |
| std | 2.54025 | 2.40738 | 0.506515 |
| min | -5.14075 | -4.29866 | -1.69449 |
| 25% | -1.97932 | -2.23894 | -0.305381 |
| 50% | 0.00854621 | 0.0123751 | 0.0778218 |
| 75% | 1.99641 | 2.16789 | 0.341274 |
| max | 5.58894 | 4.61081 | 1.44298 |
3_up.tsv
| x | y | z | |
|---|---|---|---|
| count | 38200 | 38200 | 38200 |
| mean | 3.0704e-15 | -3.60256e-15 | 1.23954e-15 |
| std | 1.69868 | 2.35114 | 0.38339 |
| min | -3.6383 | -4.07792 | -0.81214 |
| 25% | -1.53068 | -2.18586 | -0.333137 |
| 50% | -0.0697164 | 0.0415091 | 0.00216618 |
| 75% | 1.34335 | 2.10123 | 0.337469 |
| max | 3.92997 | 4.23279 | 1.12783 |
1_down.tsv
| x | y | z | |
|---|---|---|---|
| count | 38200 | 38200 | 38200 |
| mean | -1.9047e-16 | -9.52239e-15 | -2.13758e-15 |
| std | 0.272355 | 0.355092 | 0.0934928 |
| min | -2.8127 | -2.01469 | -0.868043 |
| 25% | -0.0344782 | -0.0268261 | -0.0297866 |
| 50% | -0.0105281 | -0.0028759 | 0.0181137 |
| 75% | 0.0134221 | 0.0210743 | 0.0420639 |
| max | 3.3425 | 2.05684 | 0.664769 |
2_down.tsv
| x | y | z | |
|---|---|---|---|
| count | 38200 | 38200 | 38200 |
| mean | 4.42564e-15 | 4.22792e-16 | 2.48355e-15 |
| std | 0.910304 | 0.87873 | 0.242263 |
| min | -5.11725 | -3.82467 | -1.76085 |
| 25% | -0.015856 | -0.0405354 | -0.0124892 |
| 50% | 0.00809416 | 0.007365 | 0.0354112 |
| 75% | 0.0320444 | 0.0313152 | 0.0833116 |
| max | 4.89393 | 3.9352 | 1.56822 |
3_down.tsv
| x | y | z | |
|---|---|---|---|
| count | 38200 | 38200 | 38200 |
| mean | -4.89177e-15 | -9.63697e-16 | -5.49685e-15 |
| std | 0.675041 | 0.924902 | 0.236818 |
| min | -3.52907 | -3.8726 | -3.24596 |
| 25% | -0.0323384 | -0.0405711 | -0.0126823 |
| 50% | 0.015562 | 0.00732926 | 0.0591683 |
| 75% | 0.0395122 | 0.0552297 | 0.0831185 |
| max | 3.77574 | 4.22256 | 7.81903 |
Statistical tests
- Normality test: Kolmogorov–Smirnov test
- Normality visual test: Quantile-quantile plot on chosen recording
- Stationarity test: Augmented Dickey–Fuller test
- Stationarity visual test: Autocorrelation plot
from statsmodels.tsa.stattools import adfuller
from statsmodels.api import qqplot
from scipy.stats import kstest
normality_tests = []
for name, ts in DATASET:
for x in ts.columns:
p_value = kstest(ts[x], 'norm').pvalue
test = {'name': name, 'axis': x, 'p-value': p_value, 'not-normal': p_value < 0.05}
normality_tests.append(test)
normality_tests = pd.DataFrame.from_records(normality_tests)
print(normality_tests.value_counts('not-normal'))
normality_tests.describe()
not-normal True 27 Name: count, dtype: int64
| p-value | |
|---|---|
| count | 2.700000e+01 |
| mean | 3.994415e-66 |
| std | 2.075559e-65 |
| min | 0.000000e+00 |
| 25% | 0.000000e+00 |
| 50% | 0.000000e+00 |
| 75% | 0.000000e+00 |
| max | 1.078492e-64 |
name, ts = DATASET[0]
fig, ax = plt.subplots(1, len(ts.columns), figsize=(10, 4))
for i, x in enumerate(ts.columns):
qqplot(ts[x], line='45', ax=ax[i], marker='.', alpha=0.5)
ax[i].set_title(f'Axis: {x}')
plt.tight_layout()
print(name)
plt.show()
1_still.tsv
stationarity_tests = []
for name, ts in tqdm(DATASET):
for x in ts.columns:
result = adfuller(ts[x].loc[T_WAVEFORM:T_WAVEFORM+1])
p_value = result[1]
test = {
'name': name,
'axis': x,
'statistic': result[0],
'p-value': p_value,
'stationary': p_value < 0.001
}
stationarity_tests.append(test)
stationarity_tests = pd.DataFrame.from_records(stationarity_tests)
print(stationarity_tests.value_counts('stationary'))
stationarity_tests['p-value'].describe()
0%| | 0/9 [00:00<?, ?it/s]
stationary True 26 False 1 Name: count, dtype: int64
count 2.700000e+01 mean 1.942588e-04 std 9.026014e-04 min 0.000000e+00 25% 3.192522e-29 50% 1.201432e-20 75% 4.274780e-14 max 4.677476e-03 Name: p-value, dtype: float64
name, ts = DATASET[0]
fig, ax = plt.subplots(1, len(ts.columns), figsize=(10, 4))
for i, x in enumerate(ts.columns):
ax[i].acorr(ts[x], maxlags=50)
ax[i].set_title(f'Axis: {x}')
plt.tight_layout()
print(name)
plt.show()
1_still.tsv
Time domain histogram
for name, ts in DATASET:
display(Markdown(f'**{name}**'))
axis = ts.columns
ax = ts[axis].hist(figsize=(15, 3), grid=True, bins=100, layout=(1, 3), edgecolor='black', linewidth=0.5)
plt.show()
1_still.tsv
2_still.tsv
3_still.tsv
1_up.tsv
2_up.tsv
3_up.tsv
1_down.tsv
2_down.tsv
3_down.tsv
Time domain waveform
for name, ts in DATASET:
display(Markdown(f'**{name}**'))
axis = ts.columns
ax = ts[axis].plot(figsize=(20, 8), grid=True, subplots=True)
for i, axname in enumerate(axis):
ax[i].set_xlabel('Time [s]')
ax[i].set_ylabel(f'Amplitude ({axname}) [m/s^2]')
plt.show() # plt.savefig('waveform.png')
1_still.tsv
2_still.tsv
3_still.tsv
1_up.tsv
2_up.tsv
3_up.tsv
1_down.tsv
2_down.tsv
3_down.tsv
Time domain waveform zoom detail
for name, ts in DATASET:
axis = ts.columns
display(Markdown(f'**{name}**'))
ax = (ts[axis].iloc[int(T_WAVEFORM*Fs):int(T_WAVEFORM*Fs)+Fs]
.plot(figsize=(20, 10), grid=True, subplots=True))
for i, axname in enumerate(axis):
ax[i].set_xlabel('Time [s]')
ax[i].set_ylabel(f'Amplitude ({axname}) [m/s^2]')
plt.show() # plt.savefig('waveform_zoom.png')
1_still.tsv
2_still.tsv
3_still.tsv
1_up.tsv
2_up.tsv
3_up.tsv
1_down.tsv
2_down.tsv
3_down.tsv
Time domain waveform zoom - faults side by side
fig, ax = plt.subplots(len(DATASET), 3, figsize=(12, 15), sharex=True)
for idx, df in enumerate(DATASET):
name, ts = df
columns = ts.columns
ax[idx][1].set_title(name)
ax[idx][0].set_ylabel('Amplitude [m/s^2]')
for pos, axis in enumerate(columns):
data = ts[axis].loc[T_WAVEFORM:T_WAVEFORM+0.3]
ax[idx][pos].plot(data.index, data, linewidth=1, color='darkblue')
ax[idx][pos].grid()
plt.tight_layout()
plt.show()
def spectogram(x, debug=True):
fig, ax = plt.subplots(figsize=(15, 4))
cmap = plt.get_cmap('inferno')
pxx, freqs, t, im = plt.specgram(
x, NFFT=NFFT, Fs=Fs,
detrend='mean',
mode='magnitude', scale='dB',
cmap=cmap, vmin=-60
)
fig.colorbar(im, aspect=20, pad=0.04)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Frequency [Hz]')
mafaulda.resolution_calc(Fs, NFFT)
return freqs, pxx
def window_idx(t):
return (Fs * t) // NFFT + 1
def spectrum_slice(freqs, Pxx, t):
fig, ax = plt.subplots(2, 1, figsize=(20, 8))
n = window_idx(t)
dB = 20 * np.log10(Pxx.T[n] / 0.000001)
ax[0].plot(freqs, dB) # 1 dB = 1 um/s^2
ax[0].grid(True)
ax[0].set_xlabel('Frequency [Hz]')
ax[0].set_ylabel('Amplitude [dB]')
ax[1].plot(freqs, Pxx.T[n])
ax[1].grid(True)
ax[1].set_xlabel('Frequency [Hz]')
ax[1].set_ylabel('Amplitude [m/s^2]')
return n
def get_max_frequency(freqs, Pxx, i):
max_freq = freqs[np.argmax(Pxx.T[i])]
return max_freq
def get_peaks(freqs, Pxx, i, top=5):
amplitudes = Pxx.T[i]
peaks, _ = find_peaks(amplitudes, distance=3)
fundamental = get_max_frequency(freqs, Pxx, i)
f_top = freqs[peaks[np.argsort(amplitudes[peaks])]][::-top]
y_top = np.sort(amplitudes[peaks])[::-top]
return pd.DataFrame({
'f': f_top,
'y': y_top,
'1x': f_top / fundamental
})
def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
b, a = butter(order, [lowcut, highcut], fs=fs, btype='band')
y = lfilter(b, a, data)
return y
def get_spectrograms(DATASET: List[pd.DataFrame], axis: str) -> list:
spectrograms = []
for name, ts in DATASET:
base_freq = fundamental_frequency(ts[axis], Fs)
display(Markdown(f'**{name}** *({axis.upper()} axis, Fundamental = {base_freq:.4f} Hz)*'))
freqs, Pxx = spectogram(ts[axis])
spectrograms.append((name, freqs, Pxx))
plt.show() # plt.savefig(f'x_axis_fft_{NFFT}.png')
return spectrograms
def show_spectrogram_detail(spectrograms: list, axis: str, t: float):
for name, freqs, Pxx in spectrograms:
display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
i_window = spectrum_slice(freqs, Pxx, t)
plt.show() #plt.savefig(f'x_axis_fft_{NFFT}_at_{T_SEC}s.png')
def show_mms_peaks(spectrograms: list, axis: str, t: float):
for name, freqs, Pxx in spectrograms:
display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
i_window = window_idx(t)
peaks = discovery.mms_peak_finder(Pxx.T[i_window])
fig, ax = plt.subplots(1, 1, figsize=(15, 3))
ax.grid(True)
ax.plot(freqs, Pxx.T[i_window])
ax.scatter(freqs[peaks], Pxx.T[i_window][peaks], marker='^', color='red')
ax.set_xlabel('Frequency [Hz]')
plt.show()
def show_harmonic_series(spectrograms: list, axis: str, t: float):
# https://stackoverflow.com/questions/1982770/changing-the-color-of-an-axis
for name, freqs, Pxx in spectrograms:
display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
i_window = window_idx(t)
h_series = discovery.harmonic_series_detection(freqs, Pxx.T[i_window], Fs, NFFT)
# Find best (sum of harmonics' amplitudes in the largest)
max_harmonic_amp_idx = np.argmax([
sum([h[1] for h in s]) / len(s)
for s in h_series
])
best_harmonic_series = pd.DataFrame(
h_series[max_harmonic_amp_idx],
columns=['Frequency [Hz]', 'Amplitude [m/s^2]']
)
best_harmonic_series.index += 1
display(tabulate(best_harmonic_series, headers='keys', tablefmt='html'))
# Plot found harmonic series
fig, ax = plt.subplots(1, 8, figsize=(30, 4))
for i in range(8):
s = h_series[i+1]
if i == max_harmonic_amp_idx:
ax[i].xaxis.label.set_color('red')
ax[i].plot(freqs, Pxx.T[i_window])
ax[i].scatter([x[0] for x in s], [x[1] for x in s], marker='^', color='red')
ax[i].set_xlabel('Frequency [Hz]')
plt.show()
def show_spectra_largest_amplitudes(spectrograms: list, axis: str, t: float):
for name, freqs, Pxx in spectrograms:
display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
i_window = window_idx(t)
x_fundamental = get_max_frequency(freqs, Pxx, i_window)
peaks = get_peaks(freqs, Pxx, i_window)
display(Markdown(f'- *Fundamental frequency:* {x_fundamental} Hz'))
display(tabulate(peaks.head(5), headers='keys', tablefmt='html'))
def compare_limited_specrograms(spectrograms: list, axis: str, t: float):
fig, ax = plt.subplots(len(DATASET), 1, figsize=(20, 20), sharey=True)
i = 0
for name, ts in DATASET:
signal = ts[axis].loc[t:t+NFFT/Fs].to_numpy()
pxx = np.abs(np.fft.rfft(signal) / len(signal))
freqs = np.fft.fftfreq(len(signal), d=1/Fs)[:len(pxx)]
#ilast = len(freqs[freqs < F_LIMIT])
ax[i].plot(freqs, pxx)
ax[i].grid(True)
ax[i].set_xlabel('Frequency [Hz]')
ax[i].set_ylabel('Amplitude [m/s^2]')
ax[i].set_xlim(0, F_LIMIT)
ax[i].set_ylim(0, 0.4)
ax[i].set_title(name)
i += 1
def spectrogram_energy_left_cumulative(spectrograms: list, axis: str, t: float):
fig, ax = plt.subplots(len(DATASET), 1, figsize=(20, 20), sharey=True)
i = 0
for name, ts in DATASET:
signal = ts[axis].loc[t:t+NFFT/Fs].to_numpy()
pxx = np.abs(np.fft.rfft(signal) / len(signal))
freqs = np.fft.fftfreq(len(signal), d=1/Fs)[:len(pxx)]
ax[i].plot(freqs, np.cumsum(pxx) / np.sum(pxx))
ax[i].grid(True)
ax[i].set_xlabel('Frequency [Hz]')
ax[i].set_ylabel('Cumulative energy [%]')
#ax[i].set_xlim(0, 10000)
ax[i].set_title(name)
i += 1
Compare mafaulda faults
compare_limited_specrograms(DATASET, 'x', T_SEC)
plt.tight_layout()
plt.show()
compare_limited_specrograms(DATASET, 'y', T_SEC)
plt.tight_layout()
plt.show()
compare_limited_specrograms(DATASET, 'z', T_SEC)
plt.tight_layout()
plt.show()
Compare cumulative sums
spectrogram_energy_left_cumulative(DATASET, 'x', T_SEC)
plt.tight_layout()
plt.show()
spectrogram_energy_left_cumulative(DATASET, 'y', T_SEC)
plt.tight_layout()
plt.show()
spectrogram_energy_left_cumulative(DATASET, 'z', T_SEC)
plt.tight_layout()
plt.show()
Spectrogram in X axis
x_spectra = get_spectrograms(DATASET, 'x')
1_still.tsv (X axis, Fundamental = 18.4555 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
2_still.tsv (X axis, Fundamental = 20.3534 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
3_still.tsv (X axis, Fundamental = 21.9895 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
1_up.tsv (X axis, Fundamental = 18.5209 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
2_up.tsv (X axis, Fundamental = 20.4188 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
3_up.tsv (X axis, Fundamental = 22.0550 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
1_down.tsv (X axis, Fundamental = 17.2120 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
2_down.tsv (X axis, Fundamental = 20.3534 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
3_down.tsv (X axis, Fundamental = 21.3351 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
Spectrogram detail in X axis
show_spectrogram_detail(x_spectra, 'x', T_SEC)
1_still.tsv (X axis @ 5s)
2_still.tsv (X axis @ 5s)
3_still.tsv (X axis @ 5s)
1_up.tsv (X axis @ 5s)
2_up.tsv (X axis @ 5s)
3_up.tsv (X axis @ 5s)
1_down.tsv (X axis @ 5s)
2_down.tsv (X axis @ 5s)
3_down.tsv (X axis @ 5s)
Peaks in frequency spectrum in X axis
- MMS peak finder algorithm
show_mms_peaks(x_spectra, 'x', T_SEC)
1_still.tsv (X axis @ 5s)
2_still.tsv (X axis @ 5s)
3_still.tsv (X axis @ 5s)
1_up.tsv (X axis @ 5s)
2_up.tsv (X axis @ 5s)
3_up.tsv (X axis @ 5s)
1_down.tsv (X axis @ 5s)
2_down.tsv (X axis @ 5s)
3_down.tsv (X axis @ 5s)
Harmonic series detection in X axis
# show_harmonic_series(x_spectra, 'x', T_SEC)
show_spectra_largest_amplitudes(x_spectra, 'x', T_SEC)
1_still.tsv (X axis @ 5s)
- Fundamental frequency: 19.53125 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 19.5312 | 0.546024 | 1 |
| 1 | 195.312 | 0.00964576 | 10 |
| 2 | 134.277 | 0.00462071 | 6.875 |
| 3 | 214.844 | 0.00275034 | 11 |
| 4 | 412.598 | 0.00243231 | 21.125 |
2_still.tsv (X axis @ 5s)
- Fundamental frequency: 19.53125 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 19.5312 | 1.60475 | 1 |
| 1 | 51.2695 | 0.0283099 | 2.625 |
| 2 | 195.312 | 0.010183 | 10 |
| 3 | 273.438 | 0.00566994 | 14 |
| 4 | 332.031 | 0.00470596 | 17 |
3_still.tsv (X axis @ 5s)
- Fundamental frequency: 21.97265625 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 21.9727 | 1.49894 | 1 |
| 1 | 314.941 | 0.011506 | 14.3333 |
| 2 | 329.59 | 0.00853408 | 15 |
| 3 | 97.6562 | 0.00731069 | 4.44444 |
| 4 | 9.76562 | 0.00419778 | 0.444444 |
1_up.tsv (X axis @ 5s)
- Fundamental frequency: 17.08984375 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 17.0898 | 0.468103 | 1 |
| 1 | 70.8008 | 0.0112762 | 4.14286 |
| 2 | 295.41 | 0.00670006 | 17.2857 |
| 3 | 144.043 | 0.00419564 | 8.42857 |
| 4 | 87.8906 | 0.00272608 | 5.14286 |
2_up.tsv (X axis @ 5s)
- Fundamental frequency: 19.53125 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 19.5312 | 2.06101 | 1 |
| 1 | 78.125 | 0.0223562 | 4 |
| 2 | 178.223 | 0.00997867 | 9.125 |
| 3 | 307.617 | 0.00613189 | 15.75 |
| 4 | 314.941 | 0.00486118 | 16.125 |
3_up.tsv (X axis @ 5s)
- Fundamental frequency: 21.97265625 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 21.9727 | 1.00514 | 1 |
| 1 | 195.312 | 0.0102214 | 8.88889 |
| 2 | 292.969 | 0.00891474 | 13.3333 |
| 3 | 12.207 | 0.00589523 | 0.555556 |
| 4 | 163.574 | 0.00447052 | 7.44444 |
1_down.tsv (X axis @ 5s)
- Fundamental frequency: 17.08984375 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 17.0898 | 0.0147896 | 1 |
| 1 | 53.7109 | 0.00494662 | 3.14286 |
| 2 | 205.078 | 0.00197908 | 12 |
| 3 | 83.0078 | 0.00169884 | 4.85714 |
| 4 | 173.34 | 0.00152013 | 10.1429 |
2_down.tsv (X axis @ 5s)
- Fundamental frequency: 19.53125 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 19.5312 | 0.0323693 | 1 |
| 1 | 48.8281 | 0.00282413 | 2.5 |
| 2 | 163.574 | 0.00189572 | 8.375 |
| 3 | 361.328 | 0.00152447 | 18.5 |
| 4 | 80.5664 | 0.00146388 | 4.125 |
3_down.tsv (X axis @ 5s)
- Fundamental frequency: 21.97265625 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 21.9727 | 0.0823231 | 1 |
| 1 | 65.918 | 0.006334 | 3 |
| 2 | 190.43 | 0.00178213 | 8.66667 |
| 3 | 339.355 | 0.00146555 | 15.4444 |
| 4 | 283.203 | 0.00120351 | 12.8889 |
Spectrogram in Y axis
y_spectra = get_spectrograms(DATASET, 'y')
1_still.tsv (Y axis, Fundamental = 18.4555 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
2_still.tsv (Y axis, Fundamental = 20.3534 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
3_still.tsv (Y axis, Fundamental = 21.9895 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
1_up.tsv (Y axis, Fundamental = 18.5209 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
2_up.tsv (Y axis, Fundamental = 20.4188 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
3_up.tsv (Y axis, Fundamental = 22.0550 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
1_down.tsv (Y axis, Fundamental = 17.2120 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
2_down.tsv (Y axis, Fundamental = 19.5681 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
3_down.tsv (Y axis, Fundamental = 21.3351 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
Spectrogram detail in Y axis
show_spectrogram_detail(y_spectra, 'y', T_SEC)
1_still.tsv (Y axis @ 5s)
2_still.tsv (Y axis @ 5s)
3_still.tsv (Y axis @ 5s)
1_up.tsv (Y axis @ 5s)
2_up.tsv (Y axis @ 5s)
3_up.tsv (Y axis @ 5s)
1_down.tsv (Y axis @ 5s)
2_down.tsv (Y axis @ 5s)
3_down.tsv (Y axis @ 5s)
Peaks in frequency spectrum in Y axis
show_mms_peaks(y_spectra, 'y', T_SEC)
1_still.tsv (Y axis @ 5s)
2_still.tsv (Y axis @ 5s)
3_still.tsv (Y axis @ 5s)
1_up.tsv (Y axis @ 5s)
2_up.tsv (Y axis @ 5s)
3_up.tsv (Y axis @ 5s)
1_down.tsv (Y axis @ 5s)
2_down.tsv (Y axis @ 5s)
3_down.tsv (Y axis @ 5s)
Harmonic series detection in Y axis
# show_harmonic_series(y_spectra, 'y', T_SEC)
show_spectra_largest_amplitudes(y_spectra, 'y', T_SEC)
1_still.tsv (Y axis @ 5s)
- Fundamental frequency: 19.53125 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 19.5312 | 0.823792 | 1 |
| 1 | 117.188 | 0.0125524 | 6 |
| 2 | 73.2422 | 0.00764154 | 3.75 |
| 3 | 214.844 | 0.00576193 | 11 |
| 4 | 166.016 | 0.00391798 | 8.5 |
2_still.tsv (Y axis @ 5s)
- Fundamental frequency: 19.53125 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 19.5312 | 1.63041 | 1 |
| 1 | 53.7109 | 0.0175563 | 2.75 |
| 2 | 183.105 | 0.0129831 | 9.375 |
| 3 | 292.969 | 0.00721526 | 15 |
| 4 | 266.113 | 0.00562609 | 13.625 |
3_still.tsv (Y axis @ 5s)
- Fundamental frequency: 21.97265625 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 21.9727 | 1.86651 | 1 |
| 1 | 109.863 | 0.020548 | 5 |
| 2 | 393.066 | 0.0126337 | 17.8889 |
| 3 | 361.328 | 0.009655 | 16.4444 |
| 4 | 131.836 | 0.00747923 | 6 |
1_up.tsv (Y axis @ 5s)
- Fundamental frequency: 17.08984375 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 17.0898 | 0.684795 | 1 |
| 1 | 104.98 | 0.00932373 | 6.14286 |
| 2 | 75.6836 | 0.00481567 | 4.42857 |
| 3 | 275.879 | 0.00304662 | 16.1429 |
| 4 | 305.176 | 0.00227698 | 17.8571 |
2_up.tsv (Y axis @ 5s)
- Fundamental frequency: 19.53125 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 19.5312 | 1.96031 | 1 |
| 1 | 178.223 | 0.0153813 | 9.125 |
| 2 | 185.547 | 0.00966599 | 9.5 |
| 3 | 158.691 | 0.00717196 | 8.125 |
| 4 | 244.141 | 0.00633703 | 12.5 |
3_up.tsv (Y axis @ 5s)
- Fundamental frequency: 21.97265625 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 21.9727 | 1.66467 | 1 |
| 1 | 197.754 | 0.017065 | 9 |
| 2 | 109.863 | 0.0124474 | 5 |
| 3 | 75.6836 | 0.00727467 | 3.44444 |
| 4 | 97.6562 | 0.00659581 | 4.44444 |
1_down.tsv (Y axis @ 5s)
- Fundamental frequency: 2.44140625 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 2.44141 | 0.0271932 | 1 |
| 1 | 114.746 | 0.00266754 | 47 |
| 2 | 476.074 | 0.0018169 | 195 |
| 3 | 173.34 | 0.00161734 | 71 |
| 4 | 522.461 | 0.00147466 | 214 |
2_down.tsv (Y axis @ 5s)
- Fundamental frequency: 19.53125 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 19.5312 | 0.0651357 | 1 |
| 1 | 51.2695 | 0.00371688 | 2.625 |
| 2 | 131.836 | 0.00216978 | 6.75 |
| 3 | 102.539 | 0.00176165 | 5.25 |
| 4 | 549.316 | 0.00158811 | 28.125 |
3_down.tsv (Y axis @ 5s)
- Fundamental frequency: 21.97265625 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 21.9727 | 0.0480764 | 1 |
| 1 | 85.4492 | 0.00466419 | 3.88889 |
| 2 | 244.141 | 0.0020465 | 11.1111 |
| 3 | 300.293 | 0.0017056 | 13.6667 |
| 4 | 280.762 | 0.00155137 | 12.7778 |
Spectrogram in Z axis
z_spectra = get_spectrograms(DATASET, 'z')
1_still.tsv (Z axis, Fundamental = 13.4817 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
2_still.tsv (Z axis, Fundamental = 20.3534 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
3_still.tsv (Z axis, Fundamental = 21.9895 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
1_up.tsv (Z axis, Fundamental = 0.0654 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
2_up.tsv (Z axis, Fundamental = 0.0654 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
3_up.tsv (Z axis, Fundamental = 22.0550 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
1_down.tsv (Z axis, Fundamental = 0.0654 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
2_down.tsv (Z axis, Fundamental = 0.0654 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
3_down.tsv (Z axis, Fundamental = 0.0654 Hz)
Window size: 1024 Heinsenberg box Time step: 409.6 ms Frequency step: 2.44140625 Hz
Spectrogram detail in Z axis
show_spectrogram_detail(z_spectra, 'z', T_SEC)
1_still.tsv (Z axis @ 5s)
2_still.tsv (Z axis @ 5s)
3_still.tsv (Z axis @ 5s)
1_up.tsv (Z axis @ 5s)
2_up.tsv (Z axis @ 5s)
3_up.tsv (Z axis @ 5s)
1_down.tsv (Z axis @ 5s)
2_down.tsv (Z axis @ 5s)
3_down.tsv (Z axis @ 5s)
Peaks in frequency spectrum in Z axis
show_mms_peaks(z_spectra, 'z', T_SEC)
1_still.tsv (Z axis @ 5s)
2_still.tsv (Z axis @ 5s)
3_still.tsv (Z axis @ 5s)
1_up.tsv (Z axis @ 5s)
2_up.tsv (Z axis @ 5s)
3_up.tsv (Z axis @ 5s)
1_down.tsv (Z axis @ 5s)
2_down.tsv (Z axis @ 5s)
3_down.tsv (Z axis @ 5s)
Harmonic series detection in Z axis
# show_harmonic_series(z_spectra, 'z', T_SEC)
show_spectra_largest_amplitudes(z_spectra, 'z', T_SEC)
1_still.tsv (Z axis @ 5s)
- Fundamental frequency: 14.6484375 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 14.6484 | 0.0702139 | 1 |
| 1 | 24.4141 | 0.00735616 | 1.66667 |
| 2 | 104.98 | 0.00376746 | 7.16667 |
| 3 | 163.574 | 0.00268746 | 11.1667 |
| 4 | 231.934 | 0.0020599 | 15.8333 |
2_still.tsv (Z axis @ 5s)
- Fundamental frequency: 19.53125 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 19.5312 | 0.158409 | 1 |
| 1 | 217.285 | 0.045319 | 11.125 |
| 2 | 31.7383 | 0.0148314 | 1.625 |
| 3 | 141.602 | 0.00963871 | 7.25 |
| 4 | 334.473 | 0.00469945 | 17.125 |
3_still.tsv (Z axis @ 5s)
- Fundamental frequency: 21.97265625 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 21.9727 | 0.275604 | 1 |
| 1 | 65.918 | 0.00763671 | 3 |
| 2 | 292.969 | 0.00471308 | 13.3333 |
| 3 | 205.078 | 0.00322265 | 9.33333 |
| 4 | 183.105 | 0.00215798 | 8.33333 |
1_up.tsv (Z axis @ 5s)
- Fundamental frequency: 14.6484375 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 14.6484 | 0.0538022 | 1 |
| 1 | 26.8555 | 0.00976307 | 1.83333 |
| 2 | 61.0352 | 0.00524637 | 4.16667 |
| 3 | 236.816 | 0.00400756 | 16.1667 |
| 4 | 393.066 | 0.00208041 | 26.8333 |
2_up.tsv (Z axis @ 5s)
- Fundamental frequency: 19.53125 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 19.5312 | 0.296585 | 1 |
| 1 | 117.188 | 0.0501181 | 6 |
| 2 | 178.223 | 0.0220361 | 9.125 |
| 3 | 107.422 | 0.0153412 | 5.5 |
| 4 | 205.078 | 0.0088095 | 10.5 |
3_up.tsv (Z axis @ 5s)
- Fundamental frequency: 21.97265625 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 21.9727 | 0.211261 | 1 |
| 1 | 175.781 | 0.00694372 | 8 |
| 2 | 148.926 | 0.00357859 | 6.77778 |
| 3 | 73.2422 | 0.00299298 | 3.33333 |
| 4 | 107.422 | 0.00205858 | 4.88889 |
1_down.tsv (Z axis @ 5s)
- Fundamental frequency: 2.44140625 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 2.44141 | 0.0193371 | 1 |
| 1 | 92.7734 | 0.0044802 | 38 |
| 2 | 214.844 | 0.00232566 | 88 |
| 3 | 173.34 | 0.00186004 | 71 |
| 4 | 371.094 | 0.00147204 | 152 |
2_down.tsv (Z axis @ 5s)
- Fundamental frequency: 2.44140625 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 2.44141 | 0.0535651 | 1 |
| 1 | 65.918 | 0.00347521 | 27 |
| 2 | 144.043 | 0.0028824 | 59 |
| 3 | 153.809 | 0.00226953 | 63 |
| 4 | 187.988 | 0.00185768 | 77 |
3_down.tsv (Z axis @ 5s)
- Fundamental frequency: 14.6484375 Hz
| f | y | 1x | |
|---|---|---|---|
| 0 | 14.6484 | 0.0340945 | 1 |
| 1 | 63.4766 | 0.00477927 | 4.33333 |
| 2 | 256.348 | 0.00287893 | 17.5 |
| 3 | 180.664 | 0.00190684 | 12.3333 |
| 4 | 90.332 | 0.00168929 | 6.16667 |
Histogram
axis = ['x', 'y', 'z']
for name, ts in DATASET:
display(Markdown(f'**{name}**'))
ts[axis].hist(figsize=(10, 5), grid=True, bins=50)
plt.show()
1_still.tsv
2_still.tsv
3_still.tsv
1_up.tsv
2_up.tsv
3_up.tsv
1_down.tsv
2_down.tsv
3_down.tsv
axis = ['x', 'y', 'z']
for name, ts in DATASET:
display(Markdown(f'**{name}**'))
ts[axis].boxplot(figsize=(10, 5))
plt.show()
1_still.tsv
2_still.tsv
3_still.tsv
1_up.tsv
2_up.tsv
3_up.tsv
1_down.tsv
2_down.tsv
3_down.tsv
Orbitals of all cross sections
for name, ts in DATASET:
display(Markdown(f'**{name}**'))
fig, ax = plt.subplots(1, 3, figsize=(20, 4))
for i, col in enumerate([('x', 'y'), ('x', 'z'), ('y', 'z')]):
ax[i].scatter(ts[col[0]], ts[col[1]], s=1)
ax[i].grid(True)
ax[i].set_xlabel(col[0].upper())
ax[i].set_ylabel(col[1].upper())
ax[i].grid(True)
plt.show() # plt.savefig('orbitals.png')
1_still.tsv
2_still.tsv
3_still.tsv
1_up.tsv
2_up.tsv
3_up.tsv
1_down.tsv
2_down.tsv
3_down.tsv
Orbitals of 1x harmonic frequency
x_spectra_by_name = {spec[0]: spec for spec in x_spectra}
y_spectra_by_name = {spec[0]: spec for spec in y_spectra}
z_spectra_by_name = {spec[0]: spec for spec in z_spectra}
t = 5
space = 5
for name, ts in DATASET:
display(Markdown(f'**{name}**'))
fig, ax = plt.subplots(1, 3, figsize=(20, 4))
name, freqs, Pxx = x_spectra_by_name[name]
x_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
name, freqs, Pxx = y_spectra_by_name[name]
y_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
name, freqs, Pxx = z_spectra_by_name[name]
z_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
try:
ts['x_1x'] = butter_bandpass_filter(ts['x'], x_fundamental - space, x_fundamental + space, Fs)
ts['y_1x'] = butter_bandpass_filter(ts['y'], y_fundamental - space, y_fundamental + space, Fs)
ts['z_1x'] = butter_bandpass_filter(ts['z'], z_fundamental - space, z_fundamental + space, Fs)
except ValueError:
continue
for i, col in enumerate([('x_1x', 'y_1x'), ('x_1x', 'z_1x'), ('y_1x', 'z_1x')]):
ax[i].scatter(ts[col[0]], ts[col[1]], s=1)
ax[i].grid(True)
ax[i].set_xlabel(col[0].upper())
ax[i].set_ylabel(col[1].upper())
ax[i].grid(True)
plt.show() # plt.savefig('orbitals_1x.png')
1_still.tsv
2_still.tsv
3_still.tsv
1_up.tsv
2_up.tsv
3_up.tsv
1_down.tsv
2_down.tsv
3_down.tsv
t = 5
space = 8
for name, ts in DATASET:
display(Markdown(f'**{name}**'))
name, freqs, Pxx = x_spectra_by_name[name]
x_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
name, freqs, Pxx = y_spectra_by_name[name]
y_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
name, freqs, Pxx = z_spectra_by_name[name]
z_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
try:
x = butter_bandpass_filter(ts['x'], x_fundamental - space, x_fundamental + space, Fs)
y = butter_bandpass_filter(ts['y'], y_fundamental - space, y_fundamental + space, Fs)
z = butter_bandpass_filter(ts['z'], z_fundamental - space, z_fundamental + space, Fs)
except ValueError:
continue
ax = plt.figure().add_subplot(projection='3d')
ax.scatter(x, y, z, zdir='z', s=1, color='navy')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
ax.set_zlim(-2, 2)
ax.zaxis.labelpad = -0.7
plt.show()
1_still.tsv
2_still.tsv
3_still.tsv
1_up.tsv
2_up.tsv
3_up.tsv
1_down.tsv
2_down.tsv
3_down.tsv